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ABSTRACT 

Constraints on the validity of the hierarchical gravitational instability theory and the 
evolution of biasing are presented based upon measurements of higher order clustering 
statistics in the Deeprange Survey, a catalog of ~ 710, 000 galaxies with Iab < 24 
derived from a KPNO 4m CCD imaging survey of a contiguous 4° x 4° region. We 
compute the 3-point and 4-point angular correlation functions using a direct estimation 
for the former and the counts-in-cells technique for both. The skewness S3 decreases by 
a factor of ~ 3 — 4 as galaxy magnitude increases over the range 17 < / < 22.5 (0.1 < 
z < 0.8). This decrease is consistent with a small increase of the bias with increasing 
redshift, but not by more than a factor of 2 for the highest redshifts probed. Our results 
are strongly inconsistent, at about the 3.5 — 4 cr level, with typical cosmic string models 
in which the initial perturbations follow a non-Gaussian distribution - such models 
generally predict an opposite trend in the degree of bias as a function of redshift. We 
also find that the scaling relation between the 3-point and 4-point correlation functions 
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remains approximately invariant over the above magnitude range. The simplest model 

that is consistent with these constraints is a universe in which an initially Gaussian 
perturbation spectrum evolves under the influence of gravity combined with a low level 
of bias between the matter and the galaxies that decreases slightly from z ~ 0.8 to the 
current epoch. 

Subject headings: large-scale structure, clustering, galaxy evolution, galaxy catalogs 

1. Introduction 

The evolution of the spatial distribution of galaxies is intimately related to the physical pro- 
cesses of galaxy formation, to the initial spectrum and subsequent gravitational growth of matter 
fluctuations in the early universe, and to the global geometry of space-time. Quantifying the galaxy 
distribution is, thus, fundamental to cosmology and has dominated extragalactic astronomy for the 
past two decades. The n— point correlation functions provide a statistical toolkit that can be used 
to characterize the distribution. 

The two-point correlation function is the most widely used statistic because it provides the 
most basic measure of galaxy clustering - the departure from a pure Poisson distribution. It is 
also popular because its execution is computationally straight forward. The two-point correlation 
function is defined as the joint moment of the galaxy fluctuation field, 6g, at two different positions 

6 = e = (<5s,iM, (1) 

where () means ensemble average. The two-point correlation function yields a full description of 
a Gaussian distribution only, for which all higher order connected moments are zero by definition. 
The galaxy distribution, however, exhibits non-Gaussian behavior on small scales due to non- 
linear gravitational amplification of mass fluctuations, even if they grew from an initially Gaussian 
field. On larger scales, where the density field is well represented by linear perturbation theory, 
non-Gaussian behavior may still be present if the initial perturbation spectrum was similarly non- 
Gaussian. In addition, the process of galaxy formation is likely to introduce biases between the 
clustering properties of the dark and luminous matter. In the presence of such realities, higher 
order moments arc rcqiiired to obtain a full statistical description of the galaxy distribution and to 
provide discrimination between different biasing scenarios {e.g., , Fry &; Gaztanaga 1993, Pry 1994, 
Jing 1997). 

An accurate determination of higher order clustering statistics requires a large number of 
galaxies and, to date, the most accurate measurements have been derived from wide-area angular 
surveys, such as the APM (Szapudi et al. 1995; Gaztanaga 1994) and the EDSGC (Szapudi, Meiksin, 
&; Nichol 1996), although recent redshift surveys are now becoming large enough to make interesting 
constraints (e.g., Hoyle, Szapudi, k, Baugh 1999; Szapudi et al. 2000a). These surveys are, by 
design, limited to the study of the current epoch galaxy distribution. Deep surveys add a further 
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dimension to the exploration of clustering by enabling the study of its evolution. Ideally one would 
like to have deep, wide redshift surveys available for such analyses but it's observationally infeasible 
at the current time. Projected surveys thus still provide a unique way to study the evolution of 
clustering, especially if photometric redshifts arc available. Several deep surveys have been used to 
study the evolution of low order galaxy clustering (Lilly et al. 1995; Le Fevre et al. 1995; Neuschaffer 
&; Windhorst 1995; Campos et al. 1995; Connolly et al. 1996; Lidman Sz Peterson 1996; Woods 
&; Fahlman 1997; Connolly, Szalay, &; Brunner 1998) but most suffer from insufficient contiguous 
area and are therefore barely large enough to measure the two-point correlation function. The 
Deeprange survey (Postman et al. 1998; hereafter paper I) was designed to study the evolution of 
clustering out z ^ 1. The resulting catalog contains ~ 710, 000 galaxies with Iab < 24 derived 
from a KPNO 4m CCD imaging survey of a contiguous 4° x 4° region. The photometric calibration 
of the catalog is precise enough over the entire survey to limit zeropoint drifts to ^ 0.04 mag that 
translates to a systematic error in the angular two-point correlation function, uj{B), of ^ 0.003 on a 
4° scale and proportionally less on smaller scales (sec paper I for details). Accurate measurements 
of the /-band number counts over the range 12 < Iab < 24 and the two-point angular correlation 
function up to degree scales arc presented in paper I. 

The size and quality of the Deeprange catalog are sufficient to enable reliable estimation of the 
3-point and 4-point angular correlation functions as well. This paper presents the first attempt 
to constrain higher order correlation functions down to flux limits of Iab = 23, corresponding to 
an effective redshift of z ~ 0.75. We briefly review the astrophysics contained within the higher 
moments in section 2. Section 3 presents a description of the galaxy sample analyzed and the 
computational methods used to derive the statistics and section 4 summarizes the results. The 
implications of our calculations are presented in section 5. Technical issues concerning the fit of 
the three-point correlation are presented in Appendix A. 



2. A Brief Review of Higher Order Statistics 

We begin with a summary of the basic definitions of the higher order statistical methods used 
in this paper, highlighting their most fundamental properties. For a more in depth review, the 
reader should consult the references cited. 

The observed large-scale structures in the local universe are characterized by a high degree 
of coherence (e.g. de Lapperent, Geller, &; Huchra 1986; Shectman et al. 1996; da Costa et al. 
1998) and some features, like the "Great Wall" have undergone asymmetric gravitational collapse 
(deir Antonio, Geller, & Bothun 1996). Furthermore, the galaxy distribution traces the underlying 
dark matter in a non-linear way that may depend on time and scale. It is the non-zero higher order 
correlation functions that uniquely characterize such phenomena and allow discrimination between 
the observed galaxy distribution and a Gaussian distribution with the same variance, i. e. two-point 
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correlation function. The /^-point correlation function 

= {6162 . . . 5n) (2) 

depends on a large number of parameters (3A^ coordinates, minus 3 rotations and 3 translations). 
The number of parameters decreases if the function is integrated over part of the configuration 
space (the geometric distribution of the iV-points is often referred to as their configuration, and the 
iV-dimensional space describing it as the configuration space). 

A simple model for the AT-point correlation functions is the clustering hierarchy (e.g., Peebles 
1980) defined as 

K(N) BNk N-1 

Civ(r-i, . . . ,riv) = Qnu ^ 
where ^(r) = C2(^) = {r/ro)'^, and Qjvfe are structure constants. Their average is 

A'(jV) 

QNkBNkFNk 

where i<jvfc are the form factors associated with the shape of a cell of size unity (see Boschan, 
Szapudi & Szalay 1994 for details). For the three-point correlation function, our main concern 
here, the above form factors amount to a few percent only; if neglected then ~ 53/3 (see the 
definition of S3 below). 

If the integration domain is a particular cell of volume v, with the notation / = f /v for cell 
averaging, then the amplitude of the A/"-point correlation function can be expressed as 

Sn = ^n/^''-' = {S'')c/{SY-'- (5) 

The -Sat's, with a suitable normalization motivated by leading order perturbation theory, arc com- 
monly used to characterize the higher order, i.e. non-Gaussian, properties of the galaxy distribution 
in real surveys and N-body simulations. In the second half of the above equation the integration 
over the cell, i.e. smoothing (or filtering) is implicit. While the -Sat's do not retain all the infor- 
mation encoded in the A?^-point correlation functions, in particular their shape dependence, it is an 
extremely useful measure of clustering. It is directly related to the distribution of counts in cells (in 
the same cell v), as it is the cumulant or connected moment ^ thereof. If the shape of the cell v is 
fixed, these quantities depend only on one parameter, the size of the cell. Note the alternative no- 
tation Qn = Sn/N^~'^. These two notations differ only in their normalizations: for Qn it follows 



^Ensemble averages over the connected component have the statistical property of additivity for independent pro- 
cesses and hence receive the name "cumulants" . The alternative terminology "connected" comes from the (Feynman) 
graph representation of the statistical processes. It can be shown that cumulants correspond to graphs that have 
exactly one connected component. The mathematical definition relies on a logarithmic mapping of tlie generating 
function of ordinary moments (e.g. Szapudi & Szalay 1993) 
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from the hierarchical assumption (see later), and for Sn from perturbation theory. For Gaussian 
initial conditions with power spectrum of slope n, leading order (tree-level) perturbation theory of 
the underlying density field in an expanding universe predicts the S'tv's. For N = 3 (skewness) and 
N = A (kurtosis) the prediction is (Peebles 1980; Fry 1984; Juszkiewicz, Bouchet, &; Colombi 1993; 
Bernardeau 1994; Bernardeau 1996): 

34 

Ss =y-(n + 3) 

60712 62(n + 3) 7(n + 3)2 
=W 3 + 3 ■ 

These equations depend on scale through the local slope, n, of the initial power spectrum that, 
except for scale invariant initial conditions, varies slowly with scale in all popular cosmological 
models. According to perturbation theory and simulations, geometric corrections from J1,A are 
negligible, (e.g., Bouchet et al. 1995). The above results are valid on scales > 7h~^ Mpc. On smaller 
scales accurate measurements from iV-body simulations exist (e.g., Colombi, Szapudi, Jenkins, & 
Colberg 2000, and references therein) but the measurements are subject to a particular cosmological 
model. 

If the galaxy field, Sg, is a general non-linear function of the mass density field, S, the we can 
express this function as a Taylor series f{5) = hiS + 62(5^/2-1- .... The 6jv are the non-linear biasing 
coefficients, of which 6i = 6 is the usual bias factor connecting the two-point correlation function 
of galaxies with that of the dark matter as 

^9 = b% (7) 

For higher order correlation functions, and for the Sat's, analogous calculations relate the statistics 
of the galaxy and dark matter density fields (Kaiser 1984; Bardccn et al. 1986; Grinstcin & Wise 
1986; Matarrese, Lucchin, & Bonometto 1986; Szalay 1988; Szapudi 1994; Matsubara 1995; Fry 
1996; Szapudi 1999). For example, the result for third moment is (Fry &; Gaztaiiaga 1993) 

This formula is expected to hold on the same or larger scales as leading order (weakly non-linear) 
perturbation theory, i.e. on scales > lO/i^^ Mpc. On smaller scales, Szapudi, Colombi, Cole, Frenk, 
&; Hatton (2000) found the following phenomenological rule from N-body simulations 

^^'9 = ^2(N-2) ' 

where the function 6* = 6* (6) ~ b. For b > 1 the typical effect of biasing is that it decreases 
the ^at's. The above theory does not include stochastic effects, when the galaxy density field is 
a random function of the underlying dark matter field. Stochasticity typically introduces only a 
slight extra variance on the parameters of the theory and, thus, it will not be considered further in 
this paper. 
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3. Sample Definition and Analysis Methods 

The construction of the galaxy catalog is described in paper I. Briefly, the /-band survey 
consists of 256 overlapping CCD images. The field of view of each CCD is 16 arcminutes but the 
centers are spaced 15 arcminutes apart. The 1 arcminute overlap enables accurate astrometric 
and photometric calibration to be established over the entire survey. Objects were identified, 
photometered, and classified using automated software. For consistency, we analyze the higher 
order clustering properties in the same magnitude slices® used to measure the two-point correlation 
functions in paper I. 

3.1. Counts in Cells 

As our goal is to describe the higher order clustering of galaxies, we adopt statistical measures 
that are closely related to the moments of the underlying density field. The estimation of the higher 
order correlation amplitudes follows closely the method described in Szapudi, Meiksin, & Nichol 
(1996), which can be consulted for more details. Only the most relevant definitions are given next, 
together with the outline of the technique. 

The probability distribution of counts in cells, Pn{&), is the probability that an angular cell 
of dimension 9 contains N galaxies. The factorial moments of this distribution, = ^/'/v(-^)a: 
(where {N)^ = N(N — 1)--{N — k + 1) is the fc-th falling factorial of N), are indeed closely related 
to the moments of the underlying density field, {N){l + 6), through ((1 + 5)'^) = Fk/{N)^ (Szapudi 
& Szalay 1993). 

The most common method employed to relate the discrete nature of the observed galaxy distri- 
bution to the underlying continuous density field is known as infinitesimal Poisson sampling, which 
is effectively a shot noise subtraction technique. This corresponds to the assumption that, for an 
infinitesimal cell, the number of galaxies follows a Poisson distribution with the mean determined 
by the underlying field. This assumption must be approximate for galaxies, especially on small 
scales, because of possible halo interaction or overlap. Nevertheless, on the scales we are studying, 
Poisson sampling should be a good approximation. Moreover, even on scales where no underly- 
ing continuous process exists, factorial moments are the preferred way to deal with the inherent 
discreteness of galaxies (e.g., Szapudi Sz Szalay 1997). In what follows, the continuous version 
of the theory will be used for simplicity, and factorial moments are implicitly assumed wherever 
continuous moments are used in the spirit of the above. 

Next, the mean correlation function, ^ = (S^), and the amplitudes of the AT-point correlation 



Magnitude ranges are given in the Kron-Cousins /—band; the conversion to AB magnitudes is achieved by adding 
~ 0.5 mag to the Kron-Cousins values. 
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functions, Sn = {S'^)c/{S'^)^~^, are computed ^. The galaxy correlation function measured via 
counts in cells is actually a smoothed version of this function referred to as ^. If the galaxy 
correlation function is a power-law then the average correlation function, ^, is a power-law as well 
with the same slope but with a different (increased®) amplitude that is determined with Monte- 
Carlo integration. For square cells and galaxy correlations with the usual slopes (7 ~ 1.7), the 
ratio ^/^ is typically less than a factor of 2. The calculation of ^ is an alternative way to measure 
the two-point correlation function and it is used in the computation of the normalization of the 
higher order cumulants, i.e the Sat's. 

In the sequence of calculating the S'at's from counts in cells via factorial moments, the most 
delicate step is the actual estimation of the distribution of counts in cells in the survey. We begin by 
transforming celestial coordinates into equal area Cartesian coordinates to assure proper handling 
of the curvature of the survey boundaries on the sky. We then adopt a series of masks {e.g., 
around bright stars) to denote regions in the survey that should be excluded from analysis. These 
masks were defined during the object detection phase and help minimize the inclusion of spurious 
detections. The infinitely ovcrsampling method of (Szapudi 1997) was used to estimate Pn in 
square cells. This method completely eliminates the measurement error due to the finite number 
of sampling cells. It can be sensitive, however, to edge effects on larger scales. When all the survey 
masks are used, only small scales up to 0.16° could be studied because larger cells would nearly 
always contain an excluded region. As a sensible alternative, we eliminated all but the top 5% of the 
largest masks and repeated the analysis. This procedure extended the dynamic range of the angular 
scales probed and typically docs not significantly alter the results on smaller scales. In addition 
to shrinking the possible dynamic range available to our measurement, the large number of masks 
compromise the accuracy of the error estimate as well: the geometry of the survey (as explained in 
more detail later) could not be taken into account at the level of including the geometry, position, 
and distribution of the masks. 



3.2. Three-point correlation function 

The moments of counts in cells (i.e., the -Sjv's) are the simplest descriptors of the non- 
Gaussianity of a spatial distribution, both in terms of measurement and interpretation. However, 
there are alternative statistical measures, such as the A''-th order correlation functions, that reveal 
more information about the galaxy distribution by incorporating information about the geometry 
of higher order correlations. Indeed, the S^'s are just higher order correlation functions that have 
been integrated over a portion of the available configuration space, i.e. over the A^-points within 
a given cell. Although averaging over the configuration space suppresses noise, it also erases some 
important information that is, by contrast, kept by the A^-th order correlation functions. Specif- 



()c refers to ensemble averages over the connected component 

The amphtude is higher because f is the average value of ^ UP TO a given scale. 
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ically, the validity of the hierarchical assumption can be tested directly (e.g., the first term in 

equation Al) and ultimately the effects of bias and gravitational instability can be distinguished. 
It is therefore desirable to compute as many of the A^-th order correlation functions as the survey 
volume will allow. Three-point correlation functions have been calculated for the Zwicky, Lick, 
and Jagellonian galaxy samples by Peebles k, Groth (1975), Groth & Peebles (1977), and Peebles 
(1975), for Abell clusters by Toth, Hollosy, & Szalay (1989, hereafter THS), for the APM galaxy 
survey by Frieman & Gaztahaga (1999), and for an IRAS-sclected galaxy sample by Scoccimarro et 
al. (2000; in this case the calculation was performed in Fourier space). Only the measurements by 
Peebles and Groth were performed on scales comparable to us, the rest of the work concentrated 
on much larger angular scales. 

A larger survey is needed to directly measure the three-point correlation function, ,^3 = C,, over 
a large dynamic range in scale than that required for the measurement of skewness, 5*3. Physically, 
this can be understood: the total number of triplets available in the survey are divided into finer 
subsets in the case of ^3 than for ^3, in order to examine a larger number of parameters. A smaller 
number of triangles available for each bin of the three-point correlation function naturally causes 
larger fluctuations and, hence, larger statistical errors. 

Szapudi &; Szalay (1998, hereafter SS) proposed a set of new estimators for the AT-point cor- 
relation functions using either Monte Carlo or grid methods. We denote these series of estimators 
SS/jAT and SSvka^, respectively. They have as their primary advantage the most efficient edge cor- 
rection of any estimator developed to date. The two-point correlation function was estimated both 
with the Landy-Szalay estimator (Landy & Szalay 1993; hereafter LS), {DD — 2DR+RR)/RR, and 
the grid indicator SSh^2- We employ both methods primarily to verify that the SSw2 results agree 
with those derived from the well-tested and widely used LS estimator. The three-point correlation 
function was estimated by the grid indicator SS^ks- 

The estimators SS^yjv {N = 2,3) were implemented as follows. If D denotes the data counts 

in a sufficiently fine grid, and W is the characteristic function of the grid, i.e. taking values of 1 
inside the valid survey boundary, and elsewhere, the SSw2 estimator for the two-point correlation 
function is 

C02 = {DD - 2DW + WW) /WW, (10) 

where DD, DW, and WW denote the data-data, data-window, and window-window correlation 
functions. This is analogous to the LS estimator. The SSw3 estimator for the three-point correlation 
function is 

uj3 = iDDD-3DDW + 3DWW -WWW)/WWW. (11) 

The key difference between the LS and SS^2 estimators is that the latter statistic is analogous 
to Euler's method instead of Monte Carlo integration. In fact, SSrn is the direct generalization 
of the relation for the two-point function, which is LS = SSi?2. For the grid estimators SSvfat the 
random catalog R is replaced with the characteristic function of the survey estimated on a grid, 
W. The difference amounts to a slight perturbation of the radial bins. 
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The main advantage of using a grid instead of a random catalog is computational speed, which 
starts to become prohibitive for the three-point function (at least with current technology computers 
and the present algorithms). The CPU time for calculating two-point correlation functions with 
the LS estimator is on the order of hours for about n = 10^ galaxies, and it roughly scales as 
•n?. For the three-point correlation function n-times more CPU power is thus needed, i.e. it would 
take ~ 10^ hours on a fast workstation to calculate the estimator. The SS^^jv series of grid 
estimators scale again as , but here n is the number of grid points. A significant speed up was 
achieved for the SSiys estimator by storing precalculated distances of the grid-point pairs in a look 
up table; population of the look up table required ~ in? operations. This did not change the scaling 
of the actual three-point estimator, but did minimize CPU time as floating point operations are 
eliminated from the main calculation. With this method, the three-point function for all magnitude 
cuts could be calculated in about a CPU day on a 256^ grid. 

For the two-point correlation function, where it was computationally feasible, we confirmed 
that SSi^2 produced results that were consistent with the LS estimator on scales larger than a 
few grid spacings, as it should. To repeat the same test for the three-point function would have 
necessitated supercomputer resources (see discussion above). This was not deemed justifiable as 
it is clear that no qualitative differences are expected for the three-point function from the slight 
perturbation of the binning introduced by the grid estimators. 

We computed lo^ for all magnitude slices up to about 1 degree scales. The star-galaxy classifi- 
cation accuracy degrades at faint magnitudes (/ ^ 22) in ways that are detailed in paper I. In brief, 
a fraction of the faintest compact galaxies are often classified as stars. The statistical correction for 
this effect is to include some faint "stars" into the galaxy catalog based upon a selection function 
that is derived from an extrapolation of the star counts at brighter (/ < 20.5) magnitudes where 
classification is extremely reliable. The effect of the correction is to restore, at a statistical level, 
the correct galaxy-to-star ratios at faint magnitudes. While this introduces some minor contamina- 
tion of the galaxy sample by stars, the stars are presumably randomly distributed on the angular 
scales being considered. Therefore, except for the highly improbable case when a star cluster is 
included, any stellar contamination dilutes the correlations by the usual factor where fg is the 
fraction of stars in the survey. At the magnitudes where the statistical classification corrections 
are required (I > 21.5), the number of stars introduced into the galaxy catalog as a result of the 
correction is estimated to be less than 10% of the galaxy population. The effect of the dilution is, 
thus, negligible. 

An alternative approach to dealing with faint object misclassification is to only include in 
the analysis those objects initially classified as galaxies. If the misclassified galaxies have the 
same clustering properties as the correctly classified galaxies then the correlation functions will 
not change significantly (small changes may occur as a result of cosmic variance). However, if the 
compact objects likely to be misclassified as stars cluster differently, their omission could introduce 
a bias. In that case, sampling the objects classified as stars statistically according to an estimated 
misclassification ratio, as discussed above, is the remedy. 
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Lastly, one can compute the correlation functions using all detected objects regardless of their 
classification. In this case, one can again estimate the "true" galaxy correlation function by extrap- 
olating the bright stellar number counts and multiplying the correlation amplitude by the factor 
^Obj/i^Obj — ^starY where Nobj is the total number of objects in a given magnitude bin and 
Nstar is an estimate of the number of stars in the same bin. 

To determine if misclassification effects are a significant source of bias, we estimated the two- 
point correlation function using all three strategies in all slices. The four brightest slices (17 < 
/ < 21) exhibit no discrepancies between any of the methods, indicating that misclassification 
is negligible in this flux range. For the deepest slice (22 < I < 22.5), the results based on the 
statistical inclusion of faint stars according to an estimate of the misclassification rate is in excellent 
agreement with the corrected correlation function derived using all detected objects. Agreement 
with the results using only objects classified as galaxies was marginal. The 21 < / < 22 slice is 
an intermediate case. These results suggests that correction for misclassification is essential in this 
survey for I ^ 21.5. 

For the deepest slices, we have thus performed the estimation of the three-point correlation 
functions using the first statistical correction discussed above (results from this method are denoted 
with the letter "s" in Tables 1 and 3) as well as an estimation based on using only those objects 
classified initially as galaxies. The difference in the results between the two methods provides an 
upper estimate of the systematic error introduced by star-galaxy misclassification. For the two 
deepest slices (21 < / < 22, 22 < I < 22.5), the statistical procedure should yield the correct 
results. 

4. Results 

4.1. Counts in Cells 

The S3 and S4 results are shown in Figure 1. In what follows lower case letters, such as 
sn, denote the observed angular quantities and uppercase characters, such as Sn, represent the 
deprojected, spatial statistics. In each plot, open symbols show the exact estimates when all 
cells containing survey masks are excluded (hereafter referred to as the measurement type "E") 
and closed symbols show the values when only cells containing the largest 5% of the masks were 
excluded (hereafter referred to as measurement type "M"). The magnitude cuts are shown in the 
lower right hand corner of each plot. The triangular symbols refer to S3 while the rectangular ones 
to S4. The Sat data (iV = 3, 4) and their estimated uncertainties are given in Table 1. The results 
in Table 1 for I < 21 are for the "M" type measurements whereas those for / > 21 are for the 
"E" type measurements. We provide the "E" type measurements for the faintest slices because 
they provide the best fidelity and accuracy in the presence of any faint end systematic errors in the 
catalog at the expense of angular scale dynamic range (see discussion below). 
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For the scales considered, deprojection was performed using the hierarchical Limber's equation: 

SN = RnSn- (12) 

(see e.g., Szapudi, Meiksin, & Nichol 1996). Table 2 lists the R^'s for a typical choice of the 
luminosity function. The /—band luminosity function is assumed to be the same evolving Schcchter 
function as in paper I with M*{z) = —20.9 + 5logh — (3z, a = —1.1. Our fiducial parameters for 
the deprojection are (3 = 1.5, and 7 = 1.75 for the slope of the two-point spatial correlation 
function. The cosmological parameters adopted are h = 0.65, J7o = 0.3, and A = 0. The results 
are independent of the normalization of the luminosity function, and robust against variations in 
these parameters within the accepted range, or even beyond. For instance, the effect of varying all 
parameters simultaneously over the ranges /3 = 1 — 2, 7 = 1.5 — 1.75, h = 0.5 — 0.75, = 0.3 — 1 
is a ~ 10 — 15% change in i?3. A non-zero cosmological constant A was not considered; the effect 
of a reasonable non-zero A on our estimates of the higher order correlation statistics is expected 
to be bracketed by the changes in J7 tested above. Table 2 also gives i?3 for /3 = 1 and for O = 1 
to explicitly demonstrate the relatively weak dependence of the projection coefficients on these 
parameters. 

The problem with the "E" type measurement is that the large number of excluded cells severely 
limits the maximum scale for which we can derive constraints. Indeed, on scales above 10 arcminutes 
any cell would intersect a masks with almost unit probability. With "M" type measurements, we 
can increase the maximum spatial scale by a factor of 2 at the price of contaminating the estimate 
slightly. The typical effect of such a contamination is a slight decrease in the sat's, although 
it depends on the exact spatial and size distribution of the masks and, as such, is extremely 
complicated to correct for or predict. The complex distribution and approximately power law size 
distribution of the masks prevented the application of counts in cells on larger scales, still relatively 
small compared to the characteristic size of the survey. This unfortunate phenomenon can be 
understood in terms of the masks effectively chopping up the survey into a large number of smaller 
surveys. Nevertheless, the large area of the Deeprange survey is still essential for the calculation of 
higher order statistics, as otherwise the finite volume errors would render the measurement of the 
Sat's meaningless, even on the present small scales. However, the contamination decreases as the 
scale increases, since the area-ratio of the ignored mask is decreasing, and, thus, if the "M" and 
"E" type measurements agree out to a given angular scale, then the "M" measurements on larger 
scales can be considered relatively free of contamination whereas the "E" measurements may show 
edge effects. The two faintest slices, however, exhibit some disagreement between the "E" and "M" 
measurements on the smallest scales (see Figure 1) and therefore it is the "E" measurements that 
should be considered the most reliable results in these cases. 



- 12 - 



4.1.1. Error Estimation 

The errors of the measurements (columns 4, 6, 9, and 11 in Table 1) are estimated using the full 
non-linear error calculation by (Szapudi & Colombi 1996, Szapudi, Colombi, &; Bernardeau 1999), 
obtained from the FORCE (FORtran for Cosmic Errors) package. All the parameters needed for 
such a calculation were estimated from the survey self-consistcntly. These are the perimeter and 
area of the survey, the cell area, the average galaxy count and the average two-point correlation 
function over a cell, and the average two-point correlation function over the survey (estimated from 
the integral constraint fits of paper I). The higher order Sn's up to 8 were needed as well. They were 
calculated from extended perturbation theory (Colombi et al. 1997; Szapudi, Meiksin, Sz Nichol 
1996). This theory is defined by one parameter, n^ff, which is identical to n, the local slope of the 
power spectrum on weakly non-linear scales, and is purely phenomenological on small scales. It can 
be determined from 5*3 according to Equation 9, and was measured from the average 6*3 estimated 
from the survey itself. In addition, the package uses a model for the cumulant correlators, the 
connected joint moments of counts in cells. The model by Szapudi & Szalay 1993 was chosen, as 
this was shown to provide a good description of projected data (Szapudi & Szalay 1997). 

To estimate the uncertainties in the Sn's, the FORCE package uses a perturbative procedure 
that is numerically accurate for small errors. When the fractional errors become > 1, the value 
computed is meaningless and only indicates, qualitatively, that the errors are large. In those 
instances we renormalized the plotted error to 100% and the corresponding entry in Table 1 is 
given as "L". 

The accuracy of the fully non-linear error calculation method has been thoroughly checked 
(see Szapudi, Colombi, & Bernardeau 1999; Colombi et al. 1999). However, one cautionary note is 
in order: although the calculation takes into account edge effects to second order (the ratio of the 
survey area to survey perimeter) , the use of masks introduces complex geometrical constraints and 
extra edge effects. This is the main reason why the counts in cells method is limited to fairly small 
scales; the errors associated with these effects are not accurately modeled by the FORCE. 

For our two faintest bins (/ > 21), we also provide an estimate of the fractional systematic 
errors associated with uncertainties in the statistical correction for object misclassification. These 
are given as the second error value in the parentheses in Table 1. The systematic errors are just 

the standard deviations in the mean S3 and S4 computed from 10 realizations of the statistical 
misclassification correction described in §3.2. The fractional systematic errors are shown by the 
right shifted error bars in Figure 1. Systematic errors for the brighter bins are negligible ( ^ 2%). 

4.2. Three-point correlations 

Figure 2 shows our measurements of the angular three-point correlation function, z{a, b, c), as a 
function of the hierarchical term u;{a)(jj{b)+(jj{b)uj{c)+uj{c)u{a), where u is the two-point correlation 
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function. The magnitude ranges are shown in the lower right hand corner of each subplot. The 
three-point correlation function was estimated in logarithmic bins for each side of the corresponding 
triangle. The figure makes no attempt to display shape information, although it is available from 
our estimator. The interpretation of projected shape space is extremely difficult because several 
different three-dimensional triangles project onto identical angular configurations. In the weakly 
non-linear regime, it is possible to project down the firm predictions for triangular configurations 
from a given theory and compare them with the observed angular three-point correlation function 
(e.g., Frieman & Gaztahaga 1999). For the small scales wc are considering, however, a typical shape 
would deproject to a mixture of triangular galaxy configurations that are in the highly and mildly 
non-linear regimes. Therefore, we have chosen a different, entirely phenomenological, approach 
pioneered by THS. This essentially consists of fitting the three-point correlation function with a 
class of functions motivated by a general expansion; details of the fitting procedure and error 
estimation are presented in Appendix A. 

Table 3 summarizes the results of the three-point correlation function fitting procedure. Col- 
umn 2 of this table gives the number of degrees freedom in the fit. Wc perform fits with the third 
order terms {q^^^ and q"^^) free or locked to zero - hence the two different sets of results for each mag- 
nitude bin. On scales where the correlations are small (,^2 ^ 10~^) and the relative fluctuations of 
the estimator are becoming increasingly large, the division by the hierarchical term (see Appendix 
A) becomes unstable and produces a multitude of outliers that are easily identified in Figure 2. 
Prom the measurement of ^3 from Deeprange, ^3 can be safely bracketed with 0.1 < 93 < 10, at 
least for the brighter magnitude cuts. These are fairly conservative limits and are in agreement with 
other measurements (e.g., Szapudi, Mciksin, &: Nichol 1996) as well. These limits are displayed as 
long dashes on the figures representing three-point correlation functions, and were used in the fits 
to eliminate outliers. For a few of the deepest slices the fit could be sensitive to the exact placement 
of the lower ^3 outlier cut. All the fits were, thus, reevaluated with the limits 0.01 < qs < 10 as 
well. These extended fits are given in columns 7 through 10 in Table 3. The values in columns 
3 and 7 of Table 3 are the full (unreduced) values. 

Any significant difference between the fits in Table 3 for a given magnitude interval suggests 
large systematic errors. All the fits were performed by standard computer algebra packages, and 
the parameters were found to be robust with respect to the initial value assigned to them. The 
results in Table 3 can be compared with those in Table 1 using the fact that S3 ~ 3^3. This equation 
is not exact, because the shape of the cell influences the integral performed in the definition of S3. 
The resulting form factors amount to only ~ 2 — 3%, which is the accuracy of the approximation 
(Boschan, Szapudi & Szalay 1994). Our fitting procedure cannot yield an accurate error on the 
fitted parameters, since the input errors for the were not accurately determined. Nevertheless, 
the FORCE error bars should give a good indications from the previous tables for the expected 
statistical uncertainty since the hierarchy is approximately true. Furthermore, any differences in 
the results introduced either by neglecting the statistical correction for star /galaxy misclassification 
or by altering the ^^3 fit limits (see Table 3) should give a reliable indication of the systematic errors. 
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We did not attempt to compute systematic errors for the three-point correlation function be- 
cause it would have been very time consuming computationally. We note, however, that one can 
gauge the amplitude of any systematic errors by comparing the results with and without the sta- 
tistical misclassification correction in Table 3. Furthermore, the realization of the misclassification- 
corrected galaxy catalog used in the three-point correlation function computation has S3 values 
that are very similar to the values derived from the counts-in-cells method. 

5. Discussion and Summary 

Our constraints on the higher order clustering of the brightest galaxies in our survey are in 
good agreement with previous work. This can be seen in Figure 1 for our 17 < / < 18 subsample 
where we also display the results reported by Szapudi & Gaztahaga (1998, see also Gaztaiiaga 1994) 
for the APM 17 < bj < 20 galaxy sample and those reported by Szapudi. Mciksin & Nichol (1996) 
for the EDSGC with dash-dots, and long dash-dots, respectively. Comparison of our 17 < / < 18 
subsample with the above larger surveys is justifiable because the projection coefEcients only vary 
slowly with depth. The APM and EDSGC measurements are the present state-of-the-art for shallow 
angular surveys, to be superseded only with the Sloan Digital Sky Survey. The APM and EDSGC 
measurements agree with each other on intermediate-large scales (see Szapudi &: Gaztaiiaga 1998 
for detailed comparison of counts in cells measurements of the two surveys). On small scales, 
however, the EDSGC results are higher than that of the APM. The discrepancy between the APM 
and EDSGC results can most likely be attributed to differences between the deblending procedures 
used in the construction of the two catalogs (Szapudi & Gaztahaga 1998). Our measurements, 
especially for ^3, appear to lie mostly between the two where they disagree and consistent with them, 
although perhaps somewhat lower, where they agree. We conclude that the Deeprange counts-in- 
cells measurements are in agreement with previous estimates in shallow angular surveys, despite 
the difference in wavelength and the smaller area. Constraints from the three-point correlation 
function mirror the above results as well - our measurement of Qs = 1.57 from the 17 < / < 18 
subsample is remarkably close to that of Groth and Peebles (1977), Q3 = 1.29ib0.21 compiled from 
the Zwicky, Lick, and Jagellonian samples. 

For the collection of magnitude limited cuts, ^3 and ^4 are approximately scaling. The -Sat's 
appear to decrease with depth, which suggest a small evolution with z, even though the redshift 
distribution of each slice is fairly broad. This is in accordance with theories of structure formation 
where the initial Gaussian fluctuations grow under the influence of gravity (for perturbation theory 
see e.g., Peebles 1980; Juszkiewicz, Bouchet, & Colombi 1993; Bernardeau 1992; Bernardeau 1994); 
for A?^-body simulations see e.g., Colombi et al. 1995; Baugh, Gaztanaga, &; Efstathiou 1995; Szapudi 
et al. 2000b), and where there is a small bias. 

The wide area of the Deeprange survey enables a self-contained study of the time evolution 
of the bias - there are a statistically sufficient number of low- 2; galaxies contained in the survey 
that a self-consistent constraint can be derived. Our measurement of the the bias evolution uses a 
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simple model based on the parameterization of the evolution of the two-point correlation function 
discussed in paper I. Briefly, we parameterize the redshift dependence of the correlation function 

as 

^(r,z) = (-)-^(l + z)-3-S (13) 

where 7 and tq are the slope and correlation length, and e is a phenomenological parameter describ- 
ing the evolution (e.g. Peebles 1980; Efstathiou et al. 1991; Woods & Fahlman 1997). For O = 1 
and linear evolution, the correlation function at any z can then be compared to that at the current 
epoch through the mapping (1 + ^;)^^(r/(l + z)^z). Although this mapping is strictly for f2 = 1, 
variations in the mapping due to scenarios with $7 < 1 are well within the range of uncertainties 
due to the exclusion of all details associated with specific galaxy formation processes in this simple 
model. Under the same mapping the S'at's would be invariant. Therefore, we can define the ratio 
of the mapped correlation function to that at the current epoch as 

h{zf = (1 + zf-^-'. (14) 

Since our higher order measurements probe galaxies in the highly non-linear regime, it is appropriate 
to use Equation 9 for biasing, which yields 6*3(2;) = S^{Q) /b{z)'^ . Figure 3 displays the predicted .S3, 
normalized using the measurement from the 17 < / < 18 sample, and the actual measurements at 

6 = 0.04° (solid symbols). Redshifts are assigned to the observed data by computing the median 
z for each magnitude slice as predicted by the /3 = 1.5 LF evolution model described in §4.1. The 
predicted data were computed using 7 = 1.75, and e = —1. Strictly, should be compared at the 
same comoving scale, but the fiat scaling of the S3 allows comparison at the same angular scale 
instead. This model is clearly an oversimplification of the time evolution of the bias, and yet, it 
captures the trend presented by the data remarkably well. In addition, the time evolution of the 
bias is constrained, according to Equation 9, to be less than a factor of 2 between the current epoch 
and z ~ 0.75. 

On the other hand, the above observations are in strong contrast with expectations from non- 
Gaussian scenarios. In these models, primordial non-Gaussianity (e.g., skewness) grows in linear 
theory. For Gaussian initial conditions, the growth of the skewness is a second order effect. Thus 
according to theoretical calculations (Fry & Scherrer 1994), and simulations (Colombi 1992), the 
Sm^s should have been larger in the past in non-Gaussian models compared to their Gaussian 
counterparts. As an example, 63 is expected to be a factor of 2 larger at 2; ~ 1 than at z f» for the 
typical initial conditions in cosmic string models (Stebbins 1996, private communication; Colombi 
1992). The growth of with increasing z is expected to be even more prominent. These effects 
would have been detectable in our survey, despite the dilution effects of projection and possible 
systematic errors at the survey magnitude limit. Formally, the non-Gaussian expectation is about 

7 — 8 (T from our measurement at the highest redshift. Even generously doubling our error bars 
(which would make the previous naive biasing model based on Gaussian initial conditions a perfectly 
reasonable fit to the data in Figure 3) would still exclude typical cosmic string non-Gaussian initial 
conditions at about a 3.5 — 4 cr level. The Deeprange data, thus, strongly favor Gaussian initial 
conditions. 
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While string initial conditions axe not favored for many reasons (e.g., Pen, Seljak, & Turok 
1997; Albrecht, Battye, & Robinson 1998), our expectation is that the above arguments would hold 
for a large class of non-Gaussian models. However, it is possible to invoke biasing schemes that are in 
accord with our observations but that also mask the signature of non-Gaussian initial conditions. 
For example, strong bias at early times decreases the higher order moments, counteracting the 
effects of the initial non-Gaussianity. Later the bias naturally decreases, resulting in an increase of 
the Sn's. While perhaps such a scenario is not completely unimaginable physically, e.g. by invoking 
a strong feedback during galaxy formation, it would require unnatural fine tuning in order to assure 
that the time evolution of biasing effectively cancels the time evolution of the non-Gaussian initial 
conditions. Thus, while our results cannot rule out initial non-Gaussianity with high certainty, 
the most natural explanation is that Gaussian initial conditions of the fluctuation field grew via 
gravity. Galaxies appear to trace mass quite accurately, and the small evolution of bias predicted 
by our naive model appears to describe the trends of the data fairly well. 

The three-point correlation function appears to be hierarchical down to / = 22.5. Our estimates 
of the non-hierarchical term, Q^^, are uniformly small in amplitude (see Table 3), although a cubic 
term Q^^^ cannot be excluded with high significance. While the analysis is only approximate 
in nature (the simple error model did not attempt to estimate bin-to-bin cross-correlations), the 
inclusion of a cubic term does not typically result in a significant change in the goodness of fit. 
Gravitational instability predicts Q^^^ ~ 0, thus a small cubic term is likely to mean mild bias, 
as predicted by general bias theory. In the case of a Gaussian field with a completely general 
bias function Q^^^ = is expected under fairly general conditions (Szalay 1988). Since we find 
Q^^^ <^ 1 and Q3 > 1, the galaxy density field cannot be a biased version of a Gaussian field but, 
rather, a mildly biased version of an underlying non-Gaussian field. We emphasize that the non- 
Gaussianity being referred to here is that which is induced by non-linear gravitational amplification 
and does not refer to the nature of the initial perturbation spectrum. This interpretation is subject 
to cosmic variance on Q^^^, possible systematic errors, and possible stochasticity of the bias, all of 
which could contribute to the cubic term. 

The best fitting ^3 = 3 x Q3 ~ Q^^ from the three-point correlation measurements are also 
displayed in Figure 1 as two horizontal lines: the dotted lines show result of the full fit, while 
dashed lines display the results of the restricted fit. Despite the fact that the three-point correlation 
formula extends to larger scales than the counts in cells analysis, the best fitting Q^^ appears to be 
in excellent agreement with 6*3 obtained from counts in cells, with the exception of the 19 < / < 20 
results where there is a factor of two discrepancy. The generally good agreement between the two 
methods, however, is a further indication of the insignificance of the cubic term Q^^^. In Figure 2, 
the dashed lines show the best fitting Q^^ and the dotted lines display = S3/3 as obtained 
from the counts in cells analysis. This demonstrates the agreement between the counts in cells 
analysis and the direct three-point correlation function estimation from a different perspective. 
The disagreement seen above for the 19 < / < 20 sample is less significant here and, thus, it 
must stem from the smallest scales. The points display a slight non-hierarchical curvature as well. 
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More accurate measurements from even larger surveys will be required to assess whether this is a 
significant behavior. 

In summary we have measured moments of counts in cells and the three-point correlation func- 
tion in the Deeprange survey. These constitute the deepest higher order clustering measurements 
to date. The moments measured on small scales appear to be hierarchical, and the three-point 
function, extending to larger scales continues this hierarchy. While the cubic term resulting from 
possible bias could not be excluded with high significance, the hierarchical assumption holds to a 
good approximation. This argues that gravity is the dominant process in creating galaxy corre- 
lations with bias having a detectable but minor role. Qualitatively, models with Gaussian initial 
conditions and a small amount of biasing, which increases slightly with redshift, are favored. The 
large area of the Deeprange survey allows us to study the evolution of bias over a relatively broad 
magnitude range. We find that the bias between I— band selected galaxies and the underlying 
matter distribution increases slightly with increasing redshift (up to z 0.8) but not by more than 
a factor of 2. 

In Durham, IS was supported by the PPARC rolling grant for Extragalactic Astronomy and 
Cosmology. The FORCE (FORtran for Cosmic Errors) package can be obtained from its authors, 
S. Colombi, and IS (http://www.cita.utoronto.ca/~szapudi/istvan.html). IS would like to thank 
Alex Szalay for stimulating discussions. 
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The three-point correlation function can be expanded in terms of powers of the two-point 
correlation function. This is called a Meyer clustering expansion and it is commonly used in the 
field of statistical physics as well. The second (leading) order term in the expansion corresponds 
to the hierarchical assumption, which is predicted by gravitational hierarchy (Peebles 1980). The 
higher order terms, which are determined by the nature of the galaxy-matter biasing (Szalay 1988), 
represent corrections to the hierarchy that result in non-trivial shape dependencies. 

The Meyer clustering expansion, just like any spatial expansion, can be represented by a 
Feynman-like graphical representation, shown on Figure 4. This figure also illuminates the ter- 
minology of connected components, since they correspond to connected components of the graph 
in the pictorial representation. By definition, three-point function is a connected third order mo- 
ment: it is the extra probability of a triangle above those predicted by Poisson and Gaussian (i.e. 
two-point) terms. Consequently, it makes sense to use a connected Meyer expansion for it. 

Our aim is to project a full third order connected Meyer clustering expansion of the three 
dimensional three-point function to its angular analog. Specifically, the spatial three-point corre- 
lation function can be expressed as a third order connected expansion of the two-point correlation 
function, ^, as 



This equation is directly translated to the graph on Figure 4. This expansion contains all the possi- 
ble connected terms; THS included all possible terms, including disconnected ones. Their notation 
is modified here slightly to avoid confusion with the cumulant correlator notation introduced since 



A. Fit for the Three-point Correlation Function 



6 = C(1,2,3) 



Qii(e(i)e(2) + e(2)e(3) + e(3)e(i)) + 

Q''' (^(1)^2)^(1)) + (^(1)2^(2) + sym) . 



(Al) 



THS. 



This preprint was prepared with the AAS IM^^X macros v5.0. 
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The projected three-point correlation function can be written as 

z{a,b,c)= ^ q'uji{a,b,c), (A2) 

1=11,111,21 

where a, b, and c are the angles of a triangle on the sky. The g* terms are related to the three- 
dimensional via = RiQ^. For details see THS. The projection coefficient for the hierarchical 
term Ru = R3 is identical to that of S3. Since we have found that all other terms are consistent 
with zero (i.e. the hierarchy is a good approximation) we only deproject the hierarchical term. The 
individual terms of equation Al project down to two dimensions as (THS) 

uju (ft) b, c) = uj{a)uj{b) + u{b)uj{c) + uj{c)uj{a) 
' uj{a)u){b)u;{c)\ vr^ /180 



Will (a, 6, c) - . , , , I rjf ^3 \ 

\ a + b + c J Hl-j)-^ \ TT 

U2i{a,b,c) =1 ^ ' ^ ^ +^y^] -EFT^ ' (A3) 



where ^^(7) = dx{l + x^)~'^^'^ . The second equation is only an approximation, being exact for 
7 = 2. Using the above terms, a minimum fit to the angular three-point correlation function 
was performed with the following choice 

^2 ^ ^ {z{a,b,c) - Zmod{a,b,c))'^ ^^^^ 

a,b,c 

where z^od is equation A2. The variance, from Szapudi & Szalay (1998), is 

.= = K„.(.) = ^.^. (A5) 

Here S3 = J $(a, b, c)^, a configuration integral over the definition of the angular bin, described by 
the function $(a, 6, c) = 1 when a triangle is in the bin, and zero otherwise. 

The simplification arises since S3 = S = / and SX^ ~ DDD for our choice of the character- 
istic function describing the bin ($ = $^) used in the three-point estimator (see Szapudi & Szalay 
(1998) for details). DDD is the number of data triplets in the bin. 



A.l. Error Estimation 

The variance defined in the above equation was derived for a Poisson distribution and therefore 
accounts only for discreteness effects. Two additional error contributions, edge and finite volume 
effects (Szapudi &; Colombi 1996; Szapudi, Colombi, &: Bernardeau 1999), arising from the uneven 
weights given to data points and from the fluctuations of the universe on scales larger than the 
survey size, respectively, are not accounted for. However, edge effects are expected to be small 
for these edge corrected estimators, and it can be shown that at least partial correction for finite 
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volume effects are contained in DDD. The x^'s suggest that this ansatz for the variance is within 
factor of 2 of the truth. This is the hmit of this simple model ignoring cross correlations of different 
bins. Without complicated treatment of the full correlation matrix, it would not make sense to 
improve the above formula with the inclusion of an additional phenomenological term for finite 
volume effects to tune to the number of degrees of freedom. While the results should not 
sensitively depend on the exact choice of the error model, it should be kept in mind that our errors 
for z{a, b, c) could be off by as much as a factor of 2. The most likely sense of the offset is that our 
errors may be overestimated, as judged from the x^'s- 
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Table 1. S3 and S4 as a function of scale and /—band magnitude 
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41.71 


L 


II 


0.32 


3.68 


0.63 


12.30 


L 


II 


3.65 


0.84 


20.50 


0.36 


19-20 


0.01 


3.40 


0.04 


41.57 


0.49 


20-21 


2.13 


0.04 


13.48 


L 


II 


0.02 


3.65 


0.02 


41.69 


0.19 


II 


2.56 


0.02 


16.98 


0.48 


II 


0.04 


4.54 


0.02 


53.93 


0.19 


II 


2.82 


0.03 


17.91 


0.48 


II 


0.08 


4.89 


0.05 


68.70 


0.36 


II 


2.29 


0.09 


12.15 


0.44 


II 


0.16 


4.36 


0.20 


49.25 


0.57 


II 


1.11 


0.59 




L 


II 


0.32 


4.62 


L 


2.52 


L 


II 


1.06 


L 




L 


21-22S 


0.01 


2.99 


(0.15,0.13) 


23.64 


(L,0.52) 


22-22. 5s 


3.56 


(L,0.44) 


129.8 


(L,0.26) 


II 


0.02 


2.80 


(0.05,0.05) 


25.98 


(L,0.19) 


II 


2.42 


(0.35,0.21) 


22.41 


(L,0.55) 


II 


0.04 


2.15 


(0.08,0.05) 


6.86 


(L,0.49) 


II 


2.38 


(0.31,0.10) 


6.45 


(L,0.58) 


II 


0.08 


1.73 


(0.33,0.12) 


5.84 


(L,0.57) 


II 


1.56 


(0.82,0.23) 


20.00 


(L,0.24) 


II 


0.16 


0.85 


(L,0.68) 




(L,L) 


II 


2.70 


(L,0.44) 


50.22 


(L,0.53) 



Note. — The relative error Asn/sn was obtained from the non-linear FORCE package, which uses a perturbative 
expansion in terms of the variances. Thus whenever the result is greater than 1, the error-calculation looses accuracy, 
although it is a sign that the measurements have low significance their. For such cases the entry for the table is 
"L" . Systematic error estimates are also provided for the two faintest slices and are given as the second entry in the 
parentheses. The results for slices with I < 21 are based on calculations that only exclude data within the top 5% 
largest survey masks. The two faintest slices are based on calculations in which all masked regions are excluded from 
analysis. 
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Table 2. Projection Coefficients: = S^/sn 



Imin Imax 




i?4 


i?3(/3 = 1) 


i?3(!^ = 1) 


17-18 


1.224 


1.617 


1.204 


1.232 


18-19 


1.074 


1.212 


1.189 


1.062 


19-20 


1.020 


1.051 


1.055 


1.021 


20-21 


1.049 


1.134 


1.024 


1.053 


21-22 


1.069 


1.189 


1.051 


1.071 


22-22.5 


1.076 


1.215 


1.06.5 


1.07750 



Note. — The columns R3{(3 = 1) and i?3(f2 = 1) illustrate how much R3 changes if /? = 1 (the luminosity function 
evolution parameter - see §4.1) or fl = 1 is used. 
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Table 3. fits for the general connected third order expansion of the three-point function. 







Fits with 0.1 < 93 < 10 


Fits with 0.01 < 93 


< 10 




Nf 

J 




q' 


^111 








^111 




17- 18 


44 


11.28 


1.29 


-0.26 


0.03 


11.38 


1.30 


-0.25 


0.03 


II 


42 


12.24 


1.76 






12.30 


1.76 






18 - 19 


52 


15.03 


0.86 


0.48 


0.01 


15.03 


0.86 


0.48 


0.01 


II 


50 


22.98 


1.39 






22.98 


1.39 






19- 20 


62 


81.92 


3.18 


0.29 


-0.07 


88.90 


3.24 


0.75 


-0.11 


II 


60 


85.69 


2.80 






96.45 


2.76 






20 - 21 


41 


100.31 


0.91 


-2.35 


0.16 


107.64 


1.02 


-2.11 


0.10 


II 


39 


116.11 


1.00 






126.28 


0.86 






21 - 22 


17 


10.92 


0.29 


-0.27 


0.01 


12.96 


0.15 


-0.38 


0.03 


II 


15 


11.66 


0.27 






14.06 


0.18 






22 - 22.5 


29 


28.32 


0.45 


-0.57 


0.02 


36.22 


0.47 


-0.47 


0.00 


// 


27 


3i.l5 


0.-13 






45.38 


0.32 






21 - 22s 


19 


9.95 


0.51 


-0.80 


-0.00 


13.07 


0.32 


-0.81 


0.03 


// 


17 


12.30 


0.34 






15.20 


0.24 






22 - 22.5s 


42 


25.32 


0.48 


-1.69 


0.13 


29.96 


0.39 


-1.75 


0.14 


II 


40 


33.92 


0.57 






39.56 


0.52 







Note. — Results for the deepest slices, / > 21, are shown with (mag limits followed by "s") and without statistical 
corrections for star/gala:xy misclassification (see §3.2). 
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Fig. 1. — S3 and S4 estimated from the Deeprange are shown with triangular and rectangular 
symbols, respectively. Open symbols are the results when all survey masks are used. Closed 
symbols display the results when only the top 5% largest masks are used. Dotted and dash lines 
correspond to a fit of the three-point correlation function with third order terms free or locked to 
zero, respectively. The error bars are estimated using the non-linear FORCE package. Systematic 
error estimates are also shown for the faintest (/ > 21) results as right-shifted error bars. The dash- 
dot and long dash-dot curves display the results of the APM and EDSGC surveys, respectively, in 
the upper left (17 < 7 < 18) window. 
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-6 -4 -2 -6 -4 -2 

log W^Wg + WgWg + WgW^ log W^Wg + WgWg + WgWj 

Fig. 2. — The three-point correlation function z is displayed in terms of the hierarchical term 
w{l)w{2) + w{2)w{2>) + w{2>)w{\). The long dashes show the limits that were used to exclude 
outliers in the fit for the parameters of the general third order connected expansion. Dotted 
line displays the results from the counts in cells analysis, while short dashes show the from the 

x'fit. 
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z 

Fig. 3. — Constraints on the time evolution of biasing are shown. The open symbols are the 
predicted .S3 from a simple model (see §5 for details). The solid symbols display S3 measured at 
0.04°, where the FORCE error bars are the smallest. Strietly, eaeh is at a slightly different 
scale in comoving Mpc, but because of the fiat scaling of 5*3 this is not an important effect. 
The centered error bars on all data points are those derived from the FORCE package. The right 
shifted error bars for the two faintest data points are the FORCE+systematic errors (see §4.1.1 for 
details) . 
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(g)=<^(A+A+A) + 




• • 

Fig. 4. — Graphical representation of the Meyer cluster expansion. The coordinates are denoted 
by 1,2,3 in the symbol for the three-point function, while each vertex on the right hand side of the 
"equation" contain the possible tree, and one-loop connected components of the expansion. The 
lower left corner illustrates a disconnected term as an example, which is not used in our expansion. 



